# list experiment analysis

library(foreign)
library(list)

setwd("INSERT FILE PATH HERE")
dat <- read.csv("list.csv")

# test for design effect
design.effect.test <- ict.test(dat$list, dat$treat2, J = 3, alpha = .1, gms = TRUE, pi.table = TRUE)
print(design.effect.test)

#########################
### For Fitted Values ###
#########################

dat.educ1 <- dat.educ2 <- dat.educ3 <- dat.educ4 <- dat
dat.skintone1 <- dat.skintone2 <- dat.skintone3 <- dat

dat.educ1[, "educ2"] <- 1
dat.educ2[, "educ2"] <- 2
dat.educ3[, "educ2"] <- 3
dat.educ4[, "educ2"] <- 4

dat.skintone1[, "skintone2"] <- 1
dat.skintone2[, "skintone2"] <- 2
dat.skintone3[, "skintone2"] <- 3

###########################
### Difference-in-means ###
###########################

diff.in.means <- ictreg(list ~ 1, data = dat, treat="treat2", J=3, method = "lm")
diff.in.means

#############################
### ML Model, Constrained ###
#############################

# Maximum Likelihood models (constrained)
ml.results <- ictreg(list ~ age + female + educ2 + skintone2, data = dat, treat = "treat2", 
                     J=3, method = "ml", constrained = TRUE)
summary(ml.results)

ml.predict <- predict(ml.results, avg = TRUE, interval = c("confidence"), level = .9)
ml.predict

avg.pred.educ1.ml <- predict(ml.results, newdata = dat.educ1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ2.ml <- predict(ml.results, newdata = dat.educ2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ3.ml <- predict(ml.results, newdata = dat.educ3, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ4.ml <- predict(ml.results, newdata = dat.educ4, se.fit = TRUE, avg = TRUE, interval = "confidence")

avg.pred.skintone1.ml <- predict(ml.results, newdata = dat.skintone1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone2.ml <- predict(ml.results, newdata = dat.skintone2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone3.ml <- predict(ml.results, newdata = dat.skintone3, se.fit = TRUE, avg = TRUE, interval = "confidence")

###############################
### ML Model, unconstrained ###
###############################

ml.results.uc <- ictreg(list ~ age + female + educ2 + skintone2, data = dat, treat = "treat2", 
                        J=3, method = "ml", constrained = FALSE, overdispersed = FALSE)
summary(ml.results.uc)

ml.predict.uc <- predict(ml.results.uc, avg = TRUE, interval = c("confidence"), level = .9)
ml.predict.uc

avg.pred.educ1.ml.uc <- predict(ml.results.uc, newdata = dat.educ1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ2.ml.uc <- predict(ml.results.uc, newdata = dat.educ2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ3.ml.uc <- predict(ml.results.uc, newdata = dat.educ3, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ4.ml.uc <- predict(ml.results.uc, newdata = dat.educ4, se.fit = TRUE, avg = TRUE, interval = "confidence")

avg.pred.skintone1.ml.uc <- predict(ml.results.uc, newdata = dat.skintone1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone2.ml.uc <- predict(ml.results.uc, newdata = dat.skintone2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone3.ml.uc <- predict(ml.results.uc, newdata = dat.skintone3, se.fit = TRUE, avg = TRUE, interval = "confidence")


####################
### Linear Model ###
####################

lm.results <- ictreg(list ~ age + female + educ2 + skintone2, data = dat, treat = "treat2", J=3, method = "lm")
summary(lm.results)

lm.predict <- predict(lm.results, avg = TRUE, interval = c("confidence"), level = .9)
lm.predict

avg.pred.educ1.lm <- predict(lm.results, newdata = dat.educ1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ2.lm <- predict(lm.results, newdata = dat.educ2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ3.lm <- predict(lm.results, newdata = dat.educ3, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ4.lm <- predict(lm.results, newdata = dat.educ4, se.fit = TRUE, avg = TRUE, interval = "confidence")

avg.pred.skintone1.lm <- predict(lm.results, newdata = dat.skintone1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone2.lm <- predict(lm.results, newdata = dat.skintone2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone3.lm <- predict(lm.results, newdata = dat.skintone3, se.fit = TRUE, avg = TRUE, interval = "confidence")

#####################################
### Nonlinear Least Squares Model ###
#####################################

nls.results <- ictreg(list ~ age + female + educ2 + skintone2, data = dat, treat = "treat2", J=3, method = "nls")
summary(nls.results)

nls.predict <- predict(nls.results, interval = c("confidence"), level = .90, avg = TRUE)
nls.predict

avg.pred.educ1.nls <- predict(nls.results, newdata = dat.educ1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ2.nls <- predict(nls.results, newdata = dat.educ2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ3.nls <- predict(nls.results, newdata = dat.educ3, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.educ4.nls <- predict(nls.results, newdata = dat.educ4, se.fit = TRUE, avg = TRUE, interval = "confidence")

avg.pred.skintone1.nls <- predict(nls.results, newdata = dat.skintone1, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone2.nls <- predict(nls.results, newdata = dat.skintone2, se.fit = TRUE, avg = TRUE, interval = "confidence")
avg.pred.skintone3.nls <- predict(nls.results, newdata = dat.skintone3, se.fit = TRUE, avg = TRUE, interval = "confidence")
